Common tansy ecological niche analyses
Ryan Briscoe Runquist
2024-02-07
#Libraries
library(raster)
library(sp)
library(rgdal)
library(terra)
library(rgeos)
library(sf)
library(RStoolbox)
library(ggplot2)
library(gbm)
library(dismo)
library(rJava)
library(rmaxent)
library(rasterVis)
library(viridis)
library(tidyverse)
library(grid)
library(gridExtra)
library(magick)
library(plotly)
library(ecospat)
library(FactoMineR)
library(factoextra)Population geography
Populations (n = 174) were chosen to span the entire distribution of Common tansy throughout the state of Minnesota. Common tansy is much more common in the northern half of the state, which is reflected in our sampling design. In the southern portion of the state, we attempted to collect plant tissue from as many reported populations as possible.
#Read in coordinates
pops <- read.csv("data/tansy_pop_info_mwsamples.csv")
xy <- pops[,5:6] #lat/lon of each of the populations in the dataset
coordinates(xy) <- ~ lon + lat
proj4string(xy) <- CRS("+proj=longlat +ellps=WGS84 +no_defs")
# plot(xy, xlim = c(-93, -92), ylim = c(43, 50), pch = 19, cex = 0.1)
# textplot(xy@coords[,1], xy@coords[,2], words = rownames(tansy_tab), cex = 0.25, new = FALSE)
pops$geo_grp <- as.factor(pops$geo_grp)
pops$ECS_fac2 <- as.factor(pops$ECS_fac2)
pops_df <- pops
coordinates(pops) <- ~ lon + lat
proj4string(pops) <- CRS("+proj=longlat +ellps=WGS84 +no_defs")
pops_df$finestructure_grp <- as.factor(pops_df$finestructure_grp)
pops_df$fs_clust_4 <- as.factor(pops_df$fs_clust_4)
#cluster 1 pops
clust1 <- which(pops_df$cs.K2.1>=0.5)
xy_1 <- xy[clust1]
#cluster 2 pops
clust2 <- which(pops_df$cs.K2.1<=0.5)
xy_2 <- xy[clust2]Sample Eco-Regions
The MN Department of Natural Resources has divided the state using the Ecological Classification System using the National Hierarchical Framework of Ecological Units (ECS website). We assigned each population to a geographical grouping as indicated by the colors in the map. Groupings cluster populations that were similar in geography and ecology.
midwest <- map_data("county", "minnesota")
midwest <- midwest[midwest$long > -100 & midwest$long < -89.5 &
midwest$lat > 41 & midwest$lat < 52,]
MNsites <- ggplot() + geom_path(data = midwest, mapping = aes(x = long, y = lat, group = group), color = "gray") +
geom_point(data = pops_df, mapping = aes(x = lon, y = lat, col = ECS_SUBSEC, group = name)) +
scale_color_manual(values = ECS_cols) +
coord_fixed(ratio = 1.5) +
theme_bw()
ggplotly(MNsites)Tansy spatially-explicit Genomic Structure
We used a conStruct analysis of populations to determine if there was spatially-explicit structure to genomic variation (GBS dataset of 3071 polymorphic loci). The analysis for K=2 again highlighted two geographically partitioned groupings. Group 1 (red) found mostly in the northeastern portion of the state, particularly St. Louis Co., along the North Shore, and into North Central MN. Group 2 has two distinct clusters - one in the Northwest portion of the state and the other in a portion of the forests of Lake Co. in the Northeast. Southern populations tend to be an equal sample of the two major groupings.
struct2 <- magick::image_read("conStruct/spatial/spK2_250K_structure.plot.chain_1.pdf")
#struct2_rot <- image_rotate(struct2, 270)
pie2 <- magick::image_read("conStruct/spatial/spK2_250K_pie.map.chain_1.pdf")
conStructK2 <- c(struct2, pie2)
image_append(image_scale(conStructK2, "x500"))
Environmental Variables
We downloaded and processed a broad range of environmental variables that may affect tansy population survival, growth, and dispersal for use in niche analyses. More detailed information about where data was obtained and how it was processed can be found in our GitHub
Niche Overlap analyses
Soil and Climate Niche
Hill-Smith analysis on mixed data types- aggregate environment
Ordination plot
#create background points for total climate
bkgrd <- spsample(MN_poly, 1000, type='stratified', iter = 1000)
projection(bkgrd) <- CRS("+proj=longlat +ellps=WGS84 +no_defs")
names <- vector(length = length(xy))
names[clust1] <- "clust1"
names[clust2] <- "clust2"
IDs <- c(names, rep("bkgrd", length(bkgrd)))
MN_pts <- as_tibble(rbind(xy@coords, bkgrd@coords))
MN_pts$ID <- IDs
land_use_class <- raster::subset(MN_env_stack, 23)
env_varbs_analyze <- stack(MN_plant_growth_env_3_stack, land_use_class)
site_envs <- raster::extract(env_varbs_analyze, MN_pts[,1:2])
site_envs_1 <- as_tibble(site_envs)
site_envs_1$IDs <- IDs
site_envs_2 <- na.omit(site_envs_1)
site_envs_2$drainage <- as.factor(site_envs_2$drainage)
site_envs_2$land_class <- as.factor(site_envs_2$land_class)
pca_clim <- dudi.pca(site_envs_2[,1:3], scannf = F, nf=2)
ecospat.plot.contrib(contrib = pca_clim$co, eigen = pca_clim$eig)
mix_env <- dudi.hillsmith(site_envs_2[,1:9], scannf = F, nf = 5)
ecospat.plot.contrib(contrib = mix_env$co[,1:2], eigen = mix_env$eig)
ecospat.plot.contrib(contrib = mix_env$co[,2:3], eigen = mix_env$eig)
ecospat.plot.contrib(contrib = mix_env$co[,c(1,3)], eigen = mix_env$eig)
Mixed CA scores
mix.scores.global <- mix_env$li
mix.scores.cl1 <- suprow(mix_env, site_envs_2[which(site_envs_2$IDs=="clust1"), 1:9])$li
mix.scores.cl2 <- suprow(mix_env, site_envs_2[which(site_envs_2$IDs=="clust2"), 1:9])$li
mix.scores.mn <- suprow(mix_env, site_envs_2[which(site_envs_2$IDs=="bkgrd"), 1:9])$li
mix.scores.tot <- as_tibble(cbind(site_envs_2$IDs, mix.scores.global[,1:2]))
colnames(mix.scores.tot) <- c("IDs", "Axis1", "Axis2")
biplot_env <- ggplot(data = mix.scores.tot, aes(x = Axis1, y = Axis2, color = IDs, shape = IDs)) +
geom_point() +
geom_point(data = subset(mix.scores.tot, IDs == "clust1"), aes(x = Axis1, y = Axis2, color = IDs, shape = IDs)) +
geom_point(data = subset(mix.scores.tot, IDs == "clust2"), aes(x = Axis1, y = Axis2, color = IDs, shape = IDs)) +
scale_shape_manual(values = c("clust1" = 19, "clust2" = 19, "bkgrd" = 21)) +
scale_color_manual(values = c("clust1" = "green", "clust2" = "red", "bkgrd" = "gray80")) +
theme_bw()
biplot_env
Testing for niche differentiation
Calculate and plot niche space
Axis 1 v. 2
mix.grid.clim.clust1 <- ecospat.grid.clim.dyn(glob=mix.scores.global[,1:2],
glob1=mix.scores.mn[,1:2],
sp=mix.scores.cl1[,1:2], R=100,
th.sp = 0,
th.env = 0,
extend.extent = c(-1,1, -1, 1))
mix.grid.clim.clust2 <- ecospat.grid.clim.dyn(glob=mix.scores.global[,1:2],
glob1=mix.scores.mn[,1:2],
sp=mix.scores.cl2[,1:2], R=100,
th.sp = 0,
th.env = 0,
extend.extent = c(-1,1, -1, 1))
ecospat.plot.niche.dyn(mix.grid.clim.clust1, mix.grid.clim.clust2, quant = 0.1, interest = 2,
title = "Niche overlap", name.axis1 = "Axis1", name.axis2 = "Axis2",
col.unf = "#882255", col.exp = "#88CCEE",
col.stab = "#FFDE47", colZ1 = "gray50", colZ2 = "gray80",
transparency = 0)
ecospat.shift.centroids(mix.scores.cl1[,1:2], mix.scores.cl2[,1:2], mix.scores.mn[,1:2], mix.scores.mn[,1:2])
{cairo_pdf(filename = "figures/enm_niche_overlap_axes_1_2.pdf", height = 6, width = 8)
ecospat.plot.niche.dyn(mix.grid.clim.clust1, mix.grid.clim.clust2, quant = 0.1, interest = 2,
title = "Niche overlap", name.axis1 = "PC1", name.axis2 = "PC2",
col.unf = "#882255", col.exp = "#88CCEE",
col.stab = "#FFDE47", colZ1 = "gray50", colZ2 = "gray80",
transparency = 0)
ecospat.shift.centroids(mix.scores.cl1[,1:2], mix.scores.cl2[,1:2], mix.scores.mn[,1:2], mix.scores.mn[,1:2])
}
dev.off()## quartz_off_screen
## 2Axis 1 v. 3
mix.grid.clim.clust1 <- ecospat.grid.clim.dyn(glob=mix.scores.global[,c(1,3)],
glob1=mix.scores.mn[,c(1,3)],
sp=mix.scores.cl1[,c(1,3)], R=100,
th.sp = 0,
th.env = 0,
extend.extent = c(-1,1, -1, 1))
mix.grid.clim.clust2 <- ecospat.grid.clim.dyn(glob=mix.scores.global[,c(1,3)],
glob1=mix.scores.mn[,c(1,3)],
sp=mix.scores.cl2[,c(1,3)], R=100,
th.sp = 0,
th.env = 0,
extend.extent = c(-1,1, -1, 1))
ecospat.plot.niche.dyn(mix.grid.clim.clust1, mix.grid.clim.clust2, quant = 0.1, interest = 2,
title = "Niche overlap", name.axis1 = "PC1", name.axis2 = "PC3",
col.unf = "#882255", col.exp = "#88CCEE",
col.stab = "#FFDE47", colZ1 = "gray50", colZ2 = "gray80",
transparency = 0)
ecospat.shift.centroids(mix.scores.cl1[,c(1,3)], mix.scores.cl2[,c(1,3)], mix.scores.mn[,c(1,3)], mix.scores.mn[,c(1,3)])
{cairo_pdf(filename = "figures/enm_niche_overlap_axes_1_3.pdf", height = 6, width = 8)
ecospat.plot.niche.dyn(mix.grid.clim.clust1, mix.grid.clim.clust2, quant = 0.1, interest = 2,
title = "Niche overlap", name.axis1 = "PC1", name.axis2 = "PC3",
col.unf = "#882255", col.exp = "#88CCEE",
col.stab = "#FFDE47", colZ1 = "gray50", colZ2 = "gray80",
transparency = 0)
ecospat.shift.centroids(mix.scores.cl1[,c(1,3)], mix.scores.cl2[,c(1,3)], mix.scores.mn[,c(1,3)], mix.scores.mn[,c(1,3)])
}
dev.off()## quartz_off_screen
## 2Axis 2 v. 3
mix.grid.clim.clust1 <- ecospat.grid.clim.dyn(glob=mix.scores.global[,2:3],
glob1=mix.scores.mn[,2:3],
sp=mix.scores.cl1[,2:3], R=100,
th.sp = 0,
th.env = 0,
extend.extent = c(-1,1, -1, 1))
mix.grid.clim.clust2 <- ecospat.grid.clim.dyn(glob=mix.scores.global[,2:3],
glob1=mix.scores.mn[,2:3],
sp=mix.scores.cl2[,2:3], R=100,
th.sp = 0,
th.env = 0,
extend.extent = c(-1,1, -1, 1))
ecospat.plot.niche.dyn(mix.grid.clim.clust1, mix.grid.clim.clust2, quant = 0.1, interest = 2,
title = "Niche overlap", name.axis1 = "PC2", name.axis2 = "PC3",
col.unf = "#882255", col.exp = "#88CCEE",
col.stab = "#FFDE47", colZ1 = "gray50", colZ2 = "gray80",
transparency = 0)
ecospat.shift.centroids(mix.scores.cl1[,2:3], mix.scores.cl2[,2:3], mix.scores.mn[,2:3], mix.scores.mn[,2:3])
{cairo_pdf(filename = "figures/enm_niche_overlap_axes_2_3.pdf", height = 6, width = 8)
ecospat.plot.niche.dyn(mix.grid.clim.clust1, mix.grid.clim.clust2, quant = 0.1, interest = 2,
title = "Niche overlap", name.axis1 = "PC2", name.axis2 = "PC3",
col.unf = "#882255", col.exp = "#88CCEE",
col.stab = "#FFDE47", colZ1 = "gray50", colZ2 = "gray80",
transparency = 0)
ecospat.shift.centroids(mix.scores.cl1[,2:3], mix.scores.cl2[,2:3], mix.scores.mn[,2:3], mix.scores.mn[,2:3])
}
dev.off()## quartz_off_screen
## 2Tests for overlap equivalency and similarity
#Overlap
D.overlap <- ecospat.niche.overlap(mix.grid.clim.clust1, mix.grid.clim.clust2, cor = TRUE)$D
paste(sprintf("The niche overlap is %f", D.overlap))## [1] "The niche overlap is 0.421473"#Niche equivalency
eq.test <- ecospat.niche.equivalency.test(mix.grid.clim.clust1, mix.grid.clim.clust2, rep = 1000,
intersection = 0.1,
overlap.alternative = "lower",
expansion.alternative = "higher",
stability.alternative = "lower",
unfilling.alternative = "lower")
ecospat.plot.overlap.test(eq.test, "D", "Equivalency")
ecospat.plot.overlap.test(eq.test, "expansion", "Equivalency")
ecospat.plot.overlap.test(eq.test, "stability", "Equivalency")
ecospat.plot.overlap.test(eq.test, "unfilling", "Equivalency")
sim.test <- ecospat.niche.similarity.test(mix.grid.clim.clust1, mix.grid.clim.clust2, rep=1000,
intersection = 0.1,
overlap.alternative = "lower",
expansion.alternative = "higher",
stability.alternative = "lower",
unfilling.alternative = "lower",
rand.type = 2)
ecospat.plot.overlap.test(sim.test, "D", "Similarity")
ecospat.plot.overlap.test(sim.test, "expansion", "Similarity")
ecospat.plot.overlap.test(sim.test, "stability", "Similarity")
ecospat.plot.overlap.test(sim.test, "unfilling", "Similarity")
Climatic Niche
PCA - climate
PCA plot
#create background points for total climate
bkgrd <- spsample(MN_poly, 1000, type='stratified', iter = 1000)
projection(bkgrd) <- CRS("+proj=longlat +ellps=WGS84 +no_defs")
IDs <- c(rep("site", length(xy)), rep("bkgrd", length(bkgrd)))
MN_pts <- as_tibble(rbind(xy@coords, bkgrd@coords))
MN_pts$ID <- IDs
site_envs <- raster::extract(MN_plant_growth_clim_3_stack, MN_pts[,1:2])
pca_env <- dudi.pca(site_envs, scannf = F, nf=2)
ecospat.plot.contrib(contrib = pca_env$co, eigen = pca_env$eig)
PCA scores
scores.global <- pca_env$li
scores.cl1 <- suprow(pca_env, site_envs[clust1, 1:3])$li
scores.cl2 <- suprow(pca_env, site_envs[clust2, 1:3])$li
scores.mn <- suprow(pca_env, site_envs[which(MN_pts$ID=="bkgrd"), 1:3])$li
names <- vector(length = length(xy))
names[clust1] <- "clust1"
names[clust2] <- "clust2"
IDs_spec <- c(names, rep("bkgrd", length(bkgrd)))
scores.tot <- as_tibble(cbind(IDs_spec, scores.global))
biplot_env <- ggplot(data = scores.tot, aes(x = Axis1, y = Axis2, color = IDs_spec, shape = IDs_spec)) +
geom_point() +
geom_point(data = subset(scores.tot, IDs_spec == "clust1"), aes(x = Axis1, y = Axis2, color = IDs_spec, shape = IDs_spec)) +
geom_point(data = subset(scores.tot, IDs_spec == "clust2"), aes(x = Axis1, y = Axis2, color = IDs_spec, shape = IDs_spec)) +
scale_shape_manual(values = c("clust1" = 19, "clust2" = 19, "bkgrd" = 21)) +
scale_color_manual(values = c("clust1" = "green", "clust2" = "red", "bkgrd" = "gray80")) +
theme_bw()
biplot_env
Testing for niche differentiation
Calculate and plot niche space
grid.clim.clust1 <- ecospat.grid.clim.dyn(glob=scores.global,
glob1=scores.mn,
sp=scores.cl1, R=100,
th.sp = 0)
grid.clim.clust2 <- ecospat.grid.clim.dyn(glob=scores.global,
glob1=scores.mn,
sp=scores.cl2, R=100,
th.sp = 0)
ecospat.plot.niche.dyn(grid.clim.clust1, grid.clim.clust2, quant = 0.25, interest = 2, title = "Niche overlap", name.axis1 = "PC1", name.axis2 = "PC2")
ecospat.shift.centroids(scores.cl1, scores.cl2, scores.mn, scores.mn)
Tests for overlap equivalency and similarity
#Overlap
D.overlap <- ecospat.niche.overlap(grid.clim.clust1, grid.clim.clust2, cor = TRUE)$D
paste(sprintf("The niche overlap is %f", D.overlap))## [1] "The niche overlap is 0.443794"#Niche equivalency
eq.test <- ecospat.niche.equivalency.test(grid.clim.clust1, grid.clim.clust2, rep = 1000,
intersection = 0.1,
overlap.alternative = "lower",
expansion.alternative = "higher",
stability.alternative = "lower",
unfilling.alternative = "lower")
ecospat.plot.overlap.test(eq.test, "D", "Equivalency")
sim.test <- ecospat.niche.similarity.test(grid.clim.clust1, grid.clim.clust2, rep=1000,
intersection = 0.1,
overlap.alternative = "lower",
expansion.alternative = "higher",
stability.alternative = "lower",
unfilling.alternative = "lower",
rand.type = 2)
ecospat.plot.overlap.test(sim.test, "D", "Similarity")
ecospat.plot.overlap.test(sim.test, "expansion", "Similarity")
ecospat.plot.overlap.test(sim.test, "stability", "Similarity")
ecospat.plot.overlap.test(sim.test, "unfilling", "Similarity")